//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
/*************************************************************************
 * template<class Mat_t, class EigVec_t, class EigVal_t>
 * class EsslSymSolver 
 *
 * solving ALL eigen value problems using essl routine
 * Mat_t is L packed
 *************************************************************************/
#ifndef BOOST_ESSL_SYMALLSOLVER_H
#define BOOST_ESSL_SYMALLSOLVER_H

#ifndef EIGENSPACEDIMENSION_H
#include "EigSolver/EigenSpaceDimension.h"
#endif

extern "C" {

  void dspev(int*, double*, double*, double*, int*, int*, double*, int*);
  void zhpev(int*, complex<double>*, double*, complex<double>*, int*, int*, double*, int*);

  void sspev(int*, float*, float*, float*, int*, int*, float*, int*);
  void chpev(int*, complex<float>*, float*, complex<float>*, int*, int*, float*, int*);

  void dsygv(int*, double*, int*, double*, int*, double*, double*,
             int*, int*, double*, int*);

  void dgegv(int*, double*, int*, double*, int*, double*, double*, double*,
            int*, int*, double*, int*);
}


template<class T>
class BoostEsslSymSolver {

  EigenSpaceDimension& States;       

public:

  typedef T value_type;
  typedef boost::numeric::ublas::upper                              upper;
  typedef typename boost::numeric::ublas::type_traits<T>::real_type real_type;
  typedef boost::numeric::ublas::triangular_matrix<value_type,upper> HamiltonianMatrix_t;
  typedef boost::numeric::ublas::matrix<value_type>          EigenStates_t;
  typedef boost::numeric::ublas::unbounded_array<value_type> WorkArray_t;
  typedef boost::numeric::ublas::unbounded_array<real_type>  RealArray_t;
  typedef boost::numeric::ublas::unbounded_array<int>        IndexArray_t;
  typedef typename WorkArray_t::size_type                    size_type;

  EigenSpaceDimension& States;       
  HamiltonianMatrix_t    H;
  EigenStates_t          Psi;  

  BoostEsslSymSolver(EigenSpaceDimension& dim): States(dim), naux(0) {
    H.resize(States.N, States.N); 
  }

  ~BoostEsslSymSolver() { }

  bool getEigvectors();
  bool getEigvalues();

  bool getEigvectors(HamiltonianMatrix_t& S);
  bool getEigvalues(HamiltonianMatrix_t& S);

  bool getEigvectors() {
    iopt = 1;
    Psi.resize(States.M, States.N);
    return solve(H.data().begin(), Psi.data().begin()); 
  }

  bool getEigvalues() {
    iopt = 0;
    return solve(H.data().begin(), H.data().begin()); 
  }

  bool getEigvectors(HamiltonianMatrix_t& S) {
    iopt = 1;
    Psi.resize(States.M, States.N);
    return solve(H.data().begin(), S.data().begin(), Psi.data().begin());
  }

  bool getEigvalues(HamiltonianMatrix_t& S) {
    iopt = 0;
    return solve(H.data().begin(), S.data().begin(), H.data().begin());
  }

  inline real_type eigval(size_type i) const { return eigVal[i];}

private:

  int iopt,naux,N;
  RealArray_t workd; // array of 3N+1
  RealArray_t eigVal;// eigen values

  bool solve(double* ap, double* z) {
    if(iopt == 0) check_workspace(States.N, 1); //only eigen values
    else          check_workspace(States.N, 2); //eigen values + eigen vectors
    dspev(&iopt, ap, eigVal.data().begin(), z, &N, &N, workd.data().begin(), &naux);
    return true;
  }

  bool solve(double* ap, double *bp, double* z) {
    if(iopt == 0) check_workspace(States.N, 1); //only eigen values
    else          check_workspace(States.N, 2); //eigen values + eigen vectors
    dsygv(&iopt, ap, &N, bp, &N, eigVal.data().begin(), z, &N, &N, workd.data().begin(), &naux);
    return true;
  }

  bool solve(complex<double>* ap, complex<double>* z) {
    if(iopt == 0) check_workspace(States.N, 3); //only eigen values
    else          check_workspace(States.N, 4); //eigen values + eigen vectors
    zhpev(&iopt, ap, eigVal.data().begin(), z, &N, &N, workd.data().begin(), &naux);
    return true;
  }

  bool solve(complex<double>* ap, complex<double>* bp, complex<double>* z) {
    cerr << "ESSL does not provide generalized Hermition solvers. " << endl;
    cerr << "Use other libraries, e.g., Lapack " << endl;
    return false;
  }


//   void solve(float* ap, float* a) {
//     sspev(&iopt, ap, eigVal.begin(), z, &N, &N, workd.begin(), &naux);
//   }

//   void solve(complex<float>* ap, complex<float>* a) {
//     cspev(&iopt, ap, eigVal.begin(), z, &N, &N, workd.begin(), &naux);
//   }

  void check_workspace(int n, int workd_size) {
    naux = n*workd_size;
    if(workd.size() < naux) {
      workd.resize(naux);
    }

    if(eigVal.size() < n) {
      eigVal.resize(n);
    }
  }

};

#endif

/***************************************************************************
 * $RCSfile: BoostEsslSymSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: BoostEsslSymSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

